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Capillary Condensation in Cylindrical Pores: Monte Carlo Study of the 
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When a fluid that undergoes a vapor to hquid transition in the bulk is confined to a long cylindrical pore, 
the phase transition is shifted (mostly due to surface effects at the walls of the pore) and rounded (due to 
finite size effects). The nature of the phase coexistence at the transition depends on the length of the pore: 
For very long pores the system is axially homogeneous at low temperatures. At the chemical potential where 
the transition takes place fluctuations occur between vapor-like and liquid-like states of the cylinder as a 
whole. At somewhat higher temperatures (but still far below bulk criticality) the system at phase coexistence 
is in an axially inhomogeneous multi-domain state, where long cylindrical liquid-like and vapor-like domains 
alternate. Using Monte Carlo simulations for the Ising/lattice gas model and the Asakura-Oosawa model of 
colloid-polymer mixtures the transition between these two different scenarios is characterized. It is shown 
that the density distribution changes gradually from a double-peak structure to a triple-peak shape, and the 
correlation length in axial direction (measuring the equilibrium domain length) becomes much smaller than 
the cylinder length. The (rounded) transition to the disordered phase of the fluid occurs when the axial 
correlation length has decreased to a value comparable to the cylinder diameter. It is also suggested that 
adsorption hysteresis vanishes when the transition from the simple domain state to the multi-domain state 
of the cylindrical pore occurs. We predict that the difference between the pore critical temperature and the 
hysteresis critical temperature should increase logarithmically with the length of the pore. 



I. INTRODUCTION 

The properties of both pure fluids and fluid mixtures 
confined to nanoporous and microporous materialsi"— 
have found a lot of interest recently, both from the 
point of view of various applications^^—, and also 
because phase transitions in confined geometry are 
a problem of fundamental importance in statistical 
thermodynamicsi"— ii^"— . Applications range from the 
technique to extract oil and gas from porous natural 
rocks, the use of artificial mesoporous materials such as 
various zeolithes as catalysts, "molecular sieves" to sep- 
arate fluids in fluid mixtures, and various microfluidic 
and nanofluidic devices^"— . While in some cases (e.g. 
vycor glassi2ii^) the random irregularity of the porous 
network structure is expected to lead to important phys- 
ical effects^, one can also study the idealized case of 
isolated long straight pores experimentally, both for pore 
widths on the scale of nanometers (e.g. filling fluids into 
carbon nanotube o— ' — ) and for pore widths on the scale 
of up to 150/irn (producing arrays of such pores in silicon 
wafers^^, e.g. for the purpose of characterization of DNA 
put into such pores2^, etc). 

Since a long time it is known that the vapor to liq- 
uid transition in pores is typically shifted relative to the 
condition where it occurs in the bulk: for lyophilic pore 
walls the condensation already occurs at a chemical po- 
tential where the vapor in the bulk is still undersatu- 
rated ( "capillary condensation" )lr^^'',21^ but for lyopho- 
bic pore walls the opposite effect is also possible ("cap- 
illary evaporation" )'^^"'^^. To characterize these phenom- 
ena quantitatively, however, one needs to understand 
the extent to which wetting (or drying, respectively) 



phenomena^^— exist in this restricted cylindrical geome- 
try (obviously, infinitely thick wetting or drying layers do 
not exist in narrow cylinders) . An experimentally impor- 
tant effect, that has also found a lot of theoretical atten- 
tion, is the "adsorption hysteresis" that obscures the true 
equilibrium behavior of capillary condensation in pores, 
at least over some range of parameter o"^'^^!^^" — . Another 
question concerns the understanding of critical phenom- 
ena when one reaches conditions where the density dif- 
ference between the vapor-like and liquid-like "phases" 
in the pore vanishes^""^. Here, one encounters a funda- 
mental problem of statistical mechanics, since the corre- 
lation length of density fluctuations can show unlimited 
growth only along one direction (the pore axis), but one 
does not at all expect any phase transition for quasi-one- 
dimensional systems with short-range forcea^^— . Never- 
theless, a lot of phase diagrams and critical points for var- 
ious fluids conflned in nanoscopic pores have been quoted 
in the experimental literature^"— i^"— and in the theo- 
retical worki"^i^iSi^i2&"H. This fluctuation-induced de- 
struction of the phase transition is also not seen in the- 
oretical work based on density-functional theor y^^i'^" (or 
related mean-field theories), and cannot be seen in com- 
puter simulations either, if one chooses pore lengths not 
much larger than the pore diameter, as is done in many 
caseai"— i^i^, or if one constrains fluctuations by other 
methods^i^. 

In the present work, we wish to contribute to the 
theoretical understanding of these problems, presenting 
computer simulations of two models, the Ising/lattice 
gas model confined in cylindrical geometry (as well as 
its two-dimensional analog, Ising strips of finite width), 
and the Asakura-Oosawa (AO) model for colloid-polymer 



2 



mixtures^, confined in cylinders with hard (infinitely re- 
pulsive) walls. A distinctive feature of our work is that 
we pay detailed attention to the dependence of various 
physical properties on the length L of the cylinder, con- 
fining attention to the (physically relevant) case L D 
throughout. 

Sec. 2 presents a selection of our numerical results 
for the Ising lattice, while Sec. 3 provides correspond- 
ing Monte Carlo data for the off-lattice AO model, and 
discusses the generic features of both models, interprets 
them in terms of phenomenological theoretical consider- 
ations, and draws some conclusions on pertinent experi- 
ments. Sec. 4 contains a brief summary of our work. 



II. THE ORDER PARAMETER DISTRIBUTION 
FUNCTION Pl.d(M) OF QUASI-ONE-DIMENSIONAL 
LATTICE GAS MODELS 

A. Ising strips in the L x D geometry for L 2> D 

Ising (lattice gas) models in quasi-one-dimensional ge- 
ometry have already been considered extensively in the 
literature (e.g<^iSi^"— ) but here we focus attention to 
an aspect which (to our knowledge) has not been stud- 
ied before, namely the relation between the correlation 
length ^ in the long direction of the strip and the dis- 
tribution Pl,d{M) of the magnetization per spin M in 
the system and the hysteresis behavior that one finds 
in Monte Carlo simulations applying the single spin flip 
algorithms^ (that realizes the kinetic Ising model with 
non-conserved magnetization--) . If one applies periodic 
boundary conditions in both x,y directions, the Hamil- 
tonian of the model simply is (we take the lattice spacing 
as our unit of length in this section) 

L D 

n^-jY^Y. + 1' J') + s{h3 + 1)] 

i=i j=i 

where we label the lattice sites by two indices in 
x,y directions, S{i,j) — ±1, J is the exchange energy, 
and H the (normalized) magnetic field. Here, we are 
interested in the limit L ^ oo for finite D. Note that 
lengths like L, D, £^ etc. are dimensionless in the fur- 
ther analysis. First, we summarize some exactly known 
results which are useful for our analysis: 

(i) The system does not develop a spontaneous mag- 
netization. Rather the spin correlation function for large 
distances x shows an exponential decay^iiii^l, for zero 
magnetic field, 

g{x) = {S{i,j)S{i + x,j))T oc exp(-x/£,D), x -^■ oc(2) 




FIG. 1. Sectors of size 180x5 of an Ising system with L = 10^, 
D = 5 at temperatures T (in units of J/ks) 1.1, 1.8, 2.2 and 
4.4 (from top to bottom). Up spins are shown in black and 
down spins are shown in gray. For T = 1.1 the sector was 
deliberately chosen such that it contains a domain wall in 
the center of the sector. The magnetization is zero for all 
snapshots. 



with the correlation length of the strip being given 
byS 



Cd-—^1o-^ E(-1)'^> (3) 

r— 1 

where (/3 = (kBT)-^), 

70 = 2/3J-Hlntanh(/3J) (4) 

and 



cosh7r = cosh(2/3J) coth(2/3J) - cos(r7r/D) . (5) 

Note that in the limit — >■ oo we simply get ^^^^ = 
—70. At low temperatures (/3J large) Eq. Q can be sim- 
ply approximated by the result (neglecting logarithmic 
corrections of order In D) 

\n^D~PDcT , (6) 

where a is the interfacial tension of the bulk two- 
dimensional Ising modelSi 

a = 2J - 13-^ ln[(l + cxp(-2/3J))/(l - exp(-2/3J))] . 

(7) 

Eq. ^ may simply be understood in terms of a de- 
scription of the Ising strip at low temperatures as a (di- 
lute) gas of domain walls oriented in the y-direction and 
separating large domains of opposite magnetization^i. 
Such a description is plausible when one looks at snap- 
shot pictures of the Ising strip (Fig. [1]). Eq. ^ simply 
follows when one asks at which length L of a quasi-one- 
dimensional system the free energy difference AF of a 
system with a domain wall (of free energy cost i^int) and 
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a system in a mono-domain configuration vanishes, tak- 
ing the entropy gain (In L) of putting the interface some- 
where, into accoun1)^i^ 



AF^Fi^,-(3-HnL 



Fr 



Da 



(8) 



However, Fig. [2] shows that in the temperature regime 
that is of interest for the present paper, 1.0 < fc^T/J < 
kBTc/J = 2/ln(y2 -h 1) « 2.269, Eq. ^ holds only 
qualitatively, but not quantitatively. The exact results 
{Eqs. ©-([S])} are also very useful to check that our 
Monte Carlo algorithm indeed provides a sufficiently ac- 
curate sampling: using the WolfS^ single cluster al- 
gorithm, systems of size L = 10^ were simulated for 
D = 5, 10 and 20. The (second moment) correlation 
length ^ in x-direction was then obtained from a sam- 
pling of the wave- vector-dependent susceptibility x(A;), 



x{k)^pLD{\M{k)\^), 

M{k) = (LD)-i ^ S{i, £) cxp{ik ■ 



(9) 



orienting fc in x-direction and using the smallest pos- 
sible value |fcmin| — fcmin = 27r/i, to obtain 



1 



2sin(fcmin/2) 



X(0)/x(A:min) - 1 



1/2 



(10) 



Eqs. ([9]), (fTO|) are known as an efficient method to mini- 
mize finite size effects on the estimation of the correlation 
length ^ when L and ^ are of comparable size^^, which 
is true for the case of the lowest temperatures studied, 
and then the direct estimation of from the spin cor- 
relation function {Eq. ([2])} becomes cumbersome. Fig. [5] 
shows that in this way it has been possible to "measure" 
the growth of the correlation length over 5 decades accu- 
rately. 

As is well known, the lack of better quantitative agree- 
ment between the exact result {Eqs. dS])-®} and the 
approximation based on Eqs. (Hl-dSl) can be attributed 
to the "capillary wave"— excitations of the inter- 
faces, which lead to an effective repulsive interaction be- 
tween neighboring domain walls^*^ leading to a correction 
to Eq. ([6]), which for large D gets replaced by^ 



(11) 



However, Eq. (jlip still fails in the vicinity of Tc where 
one rather finds^ 



(12) 



(ii) When for T < Tc the magnetic field is varied 
from positive to negative values the jump from positive 






FIG. 2. Correlation length plotted vs. T for L> = 5 (a), 
10 (b) and 20 (c). Full curve shows Eq. @ while broken 
curve shows the approximation Eq. (|6]). The error bars show 
Monte Carlo results that were extracted from systems with 
L = 10^ (cf. text). The dot highlights the correlation length 
of the strip at bulk criticalitySi, ^d{Tc) = AD/tv. Note the 
logarithmic scale of the ordinate. 
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(+Mo) to negative (— Mq) spontaneous magnetization, 
that would occur in the two-dimensional bulk, is slightly 
rounded. One &ndiM^ (M = (LD)-^ J2 ^U,^)) 



{M) = h[xo 



^^/[(2C.)- 



1/2 



}■ 

(13) 

The first term (xoo) in the curly brackets is just the 
susceptibility at phase coexistence in the bulk {D — >■ oo 
first, then iJ 0) for T < Tc. The second term describes 
the rounding of the transition: it extends over a region 
of fields where both terms in the square bracket of the 
denominator in Eq. ()13|) are of the same order- 



interfaces in y-direction at _ff = in its typical configu- 
ration, and the jump between ±Mo is controlled by the 
total volume LD of the strip, rounding occurring over 



H = ±kBT/{2MoLD) (18) 
and the susceptibility maximum being 



H ^±kBT/{2MaD^D) oi±D~^^^exp{~l3aD) . (14) 

The maximum value of the susceptibility x — 
d{M)/dH in the strip can then be readily obtained as 



Xmax = Xoo + 2^i?^D « D^'^ exp(/3ai?) . (15) 

As it should be, wc find that the region of the rounding 
{Eq. (fT4|)} times Xmax covers just the range ±Mo. 

The simple result for Xmax is easily interpreted in terms 
of the fiuctuation relation (for L — )■ oo; note that in the 
^ all relative distances occur twice) 



2M^Dj2{Sihmi + i,j)) 



e=o 



2M^D J dxexp[-x/^D] ^2M^D^d 



(16) 



Rather than taking correlations in the y-direction 

D 

exactly into account, ^ S{i,j) is approximated by 

MQDS{i, j) in the first approximation step in Eq. (|16p. 
Of course, Eq. (|T6l) is not to be used near Tc where 
Mo 0. 

Being interested in the effects due to the finite length 
L of the strip, we can replace g{x) = exp(— x/^u) by 
g{x) = exp(— a;/^£)) -I- exp[— (L — x)/S,d], to account for 
periodic boundary conditions. Thus instead of Eq. ([T5]) 
we then obtain 



Xr 



1 - exp(-L/^D) 



(17) 



Eq. (jl7l) is in fact a simple description of the crossover 
to the maximum susceptibility in the case of very short 
strips {L <^ ^d) where the system does not have any 



Xi] 



Xo 



2M^DL/kBT 



(19) 



While some aspects of these results were tested for the 
equivalent problem of a quantum Ising chain in a trans- 
verse field^i, we are not aware of a full test of these pre- 
dictions for the standard Ising model. Albano et al^^^ 
studied the finite size scaling of Ising strips in a I? x L 
geometry near the bulk critical temperature, demonstrat- 
ing scaling properties at constant aspect ratio D/L. An- 
other study^- considered capillary condensation in Ising 
strips with boundary fields, but considered the shift of 
the transition only, ignoring the rounding. 

When one studies phase transitions by Monte Carlo 
methods?^, the standard method of analysis is based on 
finite size scaling studies of the order parameter distribu- 
tion function Pl.d{M) and its moments2ii2^i2£. Working 
at if = 0, we can still use the same cluster algorithm as 
used for the Monte Carlo calculation of the correlation 
length (Fig. [5]), but now we are also interested in study- 
ing the effect of varying both L and D (Fig. [3]). We see 
that at low temperatures Pl,d{M) has the structure fa- 
miliar from studies in the standard square {L x L) or cube 
{L X L X L) geometry: there rather sharp peaks occur at 
±Mmax close to ±Mo(T) (cf. Fig. 0]). For comparison, 
the exact solution for the spontaneous magnetization of 
an infinite Ising latticc^^ is included. One can see that 
in this case finite size effects lead to slightly but system- 
atically larger values of the magnetization. 

At low temperatures, the region of Pl.d{M) in be- 
tween the peaks has a perfectly horizontal part. As is 
well known2§!-98, this fiat part is due to the existence of 
just two, non-interacting, interfaces crossing the system 
in y-direction. The free energy cost of creating two inter- 
faces (for pPiat ^ InL the entropic contribution where 
the interfaces at given M are placed, cf. Eq. can be 
neglected) is simply given by 2Jint = 2L)f7(T), and actu- 
ally the estimation of \b\Pl,d{M^^^) / Pl,d{Q)] « 2^^1nt 
is a useful method to numerically estimate a(T) ^^i^^ . 

However, all the above statements apply only when 
L <C ^d{T), and since C_d(T) decreases rapidly when T 
increases (Fig. [2]) the crossover when L and ^d{T) are 
of the same order needs to be considered. In Pl^b{M), 
this crossover shows up via a three-peak-structure: near 
M — Q & third peak grows and gains in weight as T 
is raised, and ultimately the peaks near ±Mo(T) have 
lost all their weight and just disappear in the tails of 
the central peak. In order to quantify this behavior, we 
define the weight of the middle peak as 
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FIG. 3. Distribution function Pl,d{M) plotted vs. M for a) 
D = 10, L = 80; b) D = 10, L = 480; c) L> = 5, L = 80; d) 
D = 5, L — 480. At M = from top to bottom: curves for 
decreasing temperatures as indicated. 




FIG. 4. Estimates for spontaneous magnetization extracted 
from the positions of the peaks of Pl,d{M). The continu- 
ous curve shows the exact solution for an infinite system^. 
Different choices of L and D are indicated. 



-t-m 



+ 1 



PL,D{M)dM/ / PL.D{M)dM (20) 



W 



where the minima of Pl.d{M) are denoted as ±m. 
Of course, at higher temperatures one always reaches a 
"spinodal temperature" Tsp{L,D) where Mmax and m 
merge, and then one no longer has a 3-peak structure, 
and Eq. (pOI loses its meaning: however, before this oc- 
curs W is practically indistinguishable from unity. We 
also emphasize that Tsp{L, D) depends on both L and D 
significantly, and like other "spinodals" it does not have 
any physical significance, for systems with short-range 
interactions like considered here^S.. 

Fig. O shows the variation of W with temperature for 
two choices of D and a range of values for L. We rec- 
ognize a gradual increase of W from W = (two-peak 
structure with perfectly flat variation of Pl.d{M) near 
M = 0) to 14^ = 1 (single Gaussian peak centered at 
M = 0) as T increases. However, the larger L becomes 
the more this gradual transition is depressed to lower 
temperature, and the sharper it becomes. It is interest- 
ing to correlate this transition with the fact that £,d{T) 
decreases from values where S,d(T) exceeds L to values 
where S,d(T) is much smaller than L. Thus, we have 
marked three temperatures for each curve where (T) — 
L {W is close to 0.1 there) and where £,d{T) = L/3 [w 
is close to 0.5 there, i.e. we arc in the center of this tran- 
sition region) and where £,d{T) = L/9 {W is close to 0.9 
there, i.e. the transition is essentially completed). Thus, 
we can define a transition temperature Tq(L,D) where 
at = the strip experiences a transition from a state 
where it is typically ordered (±Mo) to a state where it 
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^d{T) as {X{T) = [exp(2/3) - l]/[exp(-2^) + 1)], choos- 
ing henceforth units where J/kg = 1) 




FIG. 5. Weight W(T) of the central peak of Pl,d{M) plot- 
ted vs. temperature for D = 5 (a) and D = 10 (b). Various 
choices of L are included, as indicated. The symbols indicate 
the temperatures where ^ = L or ^ = L/3 or ^ = Z//9, respec- 
tively. Insert shows plots of To{L, D) vs. L, cf. text for the 
definition of To(L,D). 



is typically not uniformly ordered {{M) close to M = 0) 
although it is locally ordered (because the system is split 
into many domains of typical length £,d{T) ^ L). Hence 
we define Tq{D,L) implicitly via 



^d{To{L,D)) = L/3 , 



(21) 



and we define the temperature width AT of this tran- 
sition in terms of 



^d{To{L, D)-AT)/^d{To{L, D)+AT) = L/{L/9) = 9 . 

(22) 

At large enough L, where Tq{L,D) is so low that 
Eq. ^ is accurate, we can use Eqs. (O, © to rewrite 



(23) 



and hence Eqs. (pij) . (1^^ can be rearranged as, for 
L ^ oo {ks = 1, J = 1): 



X(To(L,i?)) -exp ^-ln(L/3)^ 
TQ{L,D)p^2D/[\n{L/i)] . (24) 

Similarly, in this limit the width AT becomes 

AT/To(T,T') « ln3/ln(T/9) . (25) 

Thus for T — > oo this transition temperature goes to 
zero, and the transition becomes gradually sharper and 
sharper, but the variations (Eqs. (P^ . (^5]) ) both are of 
order 1/ In T and hence very slow. The inset of Fig. [5)d 
shows that for the temperatures accessible for our study, 
Eq. is not yet accurate. 

It is possible to monitor this transition also in a more 
conventional way, recording either the temperature vari- 
ation of the second moment (M^) or the "susceptibility" , 
cf. Fig.il 



X' = PLD{{M') - {\M\)') 



(26) 



The use of Eq. ((26l) as an estimate for a "susceptibility" 
needs comment: of course, general statistical mechanics 
implies that x = {d{M)/dH)T = PLD{{M^) - (Mf ), 
so there does not appear any term involving the abso- 
lute value of the magnetization, and since for 77 = we 
also have (A/) = 0, % decreases monotonously with de- 
creasing temperature, and no maximum occurs. As long 
as L > ^d(T), (|M|) is small ((|M|) ^ for L ^ oo), 
and then x defined in Eq. (|26p differs from the cor- 
rect susceptibility by a constant factor (namely 1 — 2/7r). 
However, when Pl,d{M) for L < ^_d(T) just exhibits 
only two peaks at zbMmax, we have (|M|) « Mmax while 
still (M) = because of the symmetry of the distribu- 
tion against a sign change of M. Then x' a-s defined in 
Eq. (1^51) measures the width of the two peaks of the dis- 
tribution at ibMinax, while x ~ P^DM^^^. Thus the 
peak of x' is suitable to give information where the tran- 
sition from the multiple domain states at _ff = to single 
domain states in a finite strip occurs. 

As expected, both the peak positions of x' ^^nd the in- 
flection points of (M^) correlate nicely with the criterion 
that W — 0.5. The strong depression of this transition 
with increasing L is clearly seen. 

For T > To{L,D) the peaks of Pl,d{M) at ±M,„ax 
have disappeared, and a broad peak near M = re- 
mains. One can verify (Fig. [7]) that this peak is simply a 
Gaussian, 
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FIG. 6. Plot of (M^) (a,b) and x{c,d) vs. temperature, for a 
range of values of L, as indicated. Cases (a,c) refer to D = 5, 
cases (b,d) to D = 10. The asterisk in each curves marks the 
temperature at which W — 0.5. 
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FIG. 7. Distribution Pl,d{M) in the region where W « 1 
but T is still distinctly smaller than Tc, so a well-identifiable 
multi-domain configuration is observed. Broken curves show 
fit to Eq. (|27p . Insert compares the fitted value to the predic- 
tion, Eq. dm) 



PlMM) o:exp[^M'^/2{M^)] . (27) 

Noting that {cf. Eq. (HH)} (M^) = kBTxm^y,/LD we 
find in this region that 

{M^) = kBTxoo/LD + 2Ml{^D/L) . (28) 

The simulations are roughly compatible with this pre- 
diction (inset of Fig. [7]) . Finally, we draw attention to the 
temperature variation of the free energy barrier between 
the maximum of Pl,d{M) at Mmax and the minimum at 
m, 



AF = fcsTln[PL,D(M^ax)/PL,D(m)] 



(29) 



Fig.|8]shows a plot of AF/T vs. T for various choices of 
L and D. The temperatures where these barriers extrap- 
olate to zero would define the "spinodal temperatures" 
Tsp(D,L) already mentioned above, but this is not the 
point we want to make now: rather we emphasize that 
barriers AF w lO/csT are reached at temperatures far 
below Tc, where the local magnetization within a domain 
(Fig. 21) is still large. When AF becomes of the order of 
10A;bF or less, nucleation of domain walls becomes easy, 
when H is decreased and one wants to reverse the mag- 
netization in the system. To test this consideration, we 
have performed computations of the magnetization re- 
versal process of the Ising strips, using the single spin 
fiip Monte Carlo algorithm^^ to realize a (physically at 
least qualitatively realistic) dynamical evolution of the 
system (in terms of the Kinetic Ising model^*'). Starting 
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FIG. 8. Barrier AF against nucleation of interfaces in Ising 
strips plotted vs. temperature. Several choices of L and D are 
shown as indicated. The three rightmost curves correspond 
to D = 10, the three leftmost curves to D = 5. In each case, 
from left to right: L = 480, 240, 80. 




H 



FIG. 9. Magnetization of Ising strips for L = 480, D = 10 
plotted vs. field H at T = 1.6, 1.7, 1.8 and 2.0. Runs with 
decreasing H are shown as full curves, runs with increasing H 
as broken curves. While for T = 1.6 a strong hysteresis can 
still be observed, a detailed analysis shows that hysteresis 
disappears at T — 1.9. 



out at H = 0.05 we decrease the magnetic field in steps of 
AH = 0.001, equilibrating at each state for At = 2 mil- 
lion Monte Carlo steps per spin. Fig.|9]shows some exam- 
ples of the hysteresis loops that were recorded in this way. 
As expected, hysteresis loops become quickly narrow as 
the temperature is increased and when AF sa lOfc^T 
hysteresis essentially disappears completely. As a conse- 
quence, we see that the "hysteresis critical temperature" 
Tch, where hysteresis loops of our strips disappear, has 



nothing to do with Tc- (It also does not have anything 
to do with a finite size analog of T^, where no longer 
distinct domains of opposite magnetization in the strip 
can be distinguished, but one has more or less isotropic 
clusters of correlated spins of size S^d ~ Dl). Instead, 
it correlates rather well with Tq(L,D), the temperature 
where no longer uniformly ordered domains (over the full 
length L of the strip) are stable. 



B. Ising cylinders without surface fields 

Now we consider the analog of Eq. ^ on the simple 
cubic rather than the square lattice, but remove all lattice 
sites with x and y coordinates (when we define the z-axis 
as the axis of the cylinder) that satisfy 

As a boundary condition, we first choose the sim- 
ple free boundary condition, i.e. interactions to "miss- 
ing spins" do not occur. Of course, due to the lattice 
structure (which does not fit to a cylindrical surface) we 
have necessarily incquivalcnt sites at the surface: i.e., for 
R — 4 any cross section of the "cylinder" is not a sphere 
bounded by a circle, but rather we have 4 spins with three 
missing neighbors each (in a positive and negative x, y di- 
rections), on next nearest neighbor sites to those sites we 
have 8 spins with one missing neighbor, and then 8 spins 
with 2 missing neighbors follow. The consequence of this 
non-uniformity of the boundary condition have not been 
studied, however, since we do not consider it to be of real 
physical interest. 

If one does not apply any "surface magnetic 
fields" iJpiffiiiSi at these boundaries the single cluster 
algorithms^ can be straightforwardly implemented for 
this problem as well, and the probability distribution 
Pl,d{M) (where D = 2R is the diameter of the cylin- 
der) and its moments can be recorded, as described in 
Sec. Ill Al Fig. [TUlfT^ show that the findings indeed are 
qualitatively similar. Of course, x' d = 3 dimensions 
is defined as (for H = 0) 

x' = mm^) - {\M\r) (31) 

where N is the total number of spins belonging to the 
"cylinder" . Unlike the two-dimensional strips (with peri- 
odic boundary conditions also in the y-direction across 
the strip) now the ordering tendency is strongly sup- 
pressed already for the smallest value of L and due to the 
missing spins at the boundary the order in the "cylinder" 
is destabilized as expected. The same effect is seen when 
one studies Lx L squares oi Lx Lx L cubes oi Lx Lx D 
films with free boundaries, as is well knowni2^"-i2£. 

We also expect in this case a simple exponential decay 
of the spin correlation function in the axial direction of 
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FIG. 10. Weight of the central peak W {Eq. ((20|} of the 
order parameter distribution for three choices of the radius R 
{R = 3, 4, 5 from left to right, in different colors) and three 
choices of L in each case: L — 600, 400, 200 from left to 
right with different line styles. The arrow shows the critical 
temperature Tc^^^- of the bulk three-dimensional Ising model. 
Note that a periodic boundary condition is only applied along 
the a-axis of the "cylinder" . 
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the cylinder, analogous to Eq. but the correspond- 
ing correlation length is not known independently. In 
analogy to Eq. ([B]), we expect that ^jy at low tempera- 
tures (and L od) simply varies exponentially with the 
cross-sectional area A of the cylinder. 



oc exp(/3 Act), A = 



(32) 



where the simple relation between A and D applies for 
off-lattice models with strictly circular cross section of the 
cylinder (in the present Ising model case, A = Nc{R), 
the number of spins in a cross sectional plane for the 
considered choice of R, i.e. iVc(3) = 29, iVc(4) = 49 and 
iVc(5) = 69). Now a is the interfacial free energy per spin 
for the three-dimensional Ising model. However, using 
the same reasoning as in Eqs. (PT|) we now obtain 
that the effective transition temperature Tq{L,D) of a 
long cylinder from a multi-domain configuration to the 
single-domain configuration is given by 



kBTQ{L,D)/a{To) = A/ln{L/3) , Llarge 



(33) 



While for T — again (j{T) 2J for planar inter- 
faces, for the temperatures of interest for the present 
study Eq. p3p is not expected to be quantitatively ac- 
curate (at not so low temperatures due to boundary ef- 
fects we expect that the actual interfacial energy Fint is 
smaller than the asymptotic estimate Aa, for the small 
radii R studied here). Nevertheless, Fig. [T^ shows that 
To{L,D) exhibits a distinct decrease with increasing L, 
as expected. 
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FIG. 11. Plot of (M^) (a,b) and x'ic,d) vs. temperature, 
for a range of values of L, as indicated. Cases (a,c) refer to 
R — 3 and cases (b,d) to i? = 4. Dot on each curve shows the 
temperature for which W = 0.5. 
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FIG. 12. To{L, D) for three choices of R (as indicated) plotted 
vs. [ln(L/3)]~' 



C. Ising cylinders with surface fields 

When the Ising model is re-interpreted as a lattice 
gas, it is natural to assume that a "free surface" bound- 
ary condition is physically caused by a confining external 
wall, which then is expected to provide an external po- 
tential to the particles adjacent to the wall. In "magnetic 
language" , such a wall potential translates to a "surface 
magnetic field" 

While in the absence of Hi the spin reversal symmetry 
of the Ising model ensures that phase coexistence between 
domains of opposite spontaneous magnetization occurs 
for bulk field iJ = 0, for iJi ^0 this symmetry is broken. 
Thus, the lattice gas model in thin film geometry with 
Hi ^ has been thoroughly discussed as a model for 
capillary condensationi^i^^i^iiSi^-iSS. 

Choosing a surface field Hi = 0.75, J — 1 that acts on 
all spins in the surface layer (i.e., spins that have "miss- 
ing neighbors") we have to carry out scans where H is 
varied to locate Hcoex{T), i.e. the field where in short 
pores phase coexistence occurs. At low temperatures, 
where for the considered choices of pore length L at phase 
coexistence only two peaks occur in the distribution 
function Pl,d{M) (cf. Fig. [T3|) an accurate sampling of 
Pl,d{M) is possible combining the standard Metropolis 
algorithms^ with successive umbrella samplingii^ meth- 
ods. Note that the single cluster algorithm^^ or the re- 
lated Swendsen-Wangiii algorithm can be extended to 
include bulk and surface fields but become very ineffi- 
cient (apart from the case where both surface and bulk 
fields become extremely small-'^-SS,) and then would not 
present any advantage. 

The location of the coexistence field Hcock{T) in the 
regime where Pl,d{M) shows only two peaks can be 
based on the "equal area rule"^^^*^ as for phase transi- 
tions in the bulk. Of course, in the present case, i.e. for 
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FIG. 13. a) Plot of Pl,d{M) vs. M for a cylinder of radius 
R = 4 and length L — 200 for Hi = 0.75 and three temper- 
atures T = 3.5,3.65 and 3.75, as indicated. The chosen field 
H = -H'coGx(r) in the bulk is also shown in the figure, b) Same 
as a), but for T = 3.58, H — —0.2279 (which was a first rough 
estimate for //coox(T)) and three choices of L, as indicated. 
The final estimates for Hcoax{T) are close to H — —0.2277 
and were found by histogram extrapolation methods, requir- 
ing an equal area rule for the two outer peaks. Note that 
Pl,d{M) no longer exhibits any symmetry with respect to a 
sign change of M . 



finite pore diameter D, we still expect that this tran- 
sition never becomes sharp, irrespectively how large L 
becomes. Thus, in the limit L 00 the susceptibility x, 
which now needs to be defined by 



(34) 



becomes a delta function (the limiting behavior for first 
order transitions^iii^) only for T — > 0. In the case of L 
finite where at low T only two peaks in P^^d^M) occur, 
we expect that (M^) - (M)^ reaches at H = i?cocx(T) 
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a maximum of order unity which we denote as Cniax(2^) 
(since the spin reversal symmetry is broken due to the 
surface field, we can no longer conclude that (M) = 
at H — Hcocx{T)). Thus we conclude that x reaches a 
maximum value 



Xmax(T) = /3iVc,„ax(r), H = i/cocx(T) 



(35) 



while in the region of fields where Pl,d (M) has a single 
maximum only (at the considered low temperature where 
^d{T) ^ L) the susceptibility x will be of order unity. 
Following the reasoning of^"'^'^ we can conclude that in this 
region the rounding of the transitions is simply given by 
the condition that 
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Since N is very large throughout our study, AH is very 
small, and i?cocx(r) in this regime of low temperatures 
is rather well defined. 

Of course, the situation becomes more subtle at higher 
temperatures, where the third peak in the distribution 
Pl,d{M) appears (Fig. [T3]). As long as the weight of 
this central peak is not yet much larger than the weight 
of the two other, sharper, peaks, we simply can ignore 
this peak and still apply the equal weight rule with re- 
spect to the two outer peaks. However, the equal weight 
rule method for estimating Hcocx{T) becomes obsolete 
when the weights of the two outer peaks become rela- 
tively small (and ultimately the two outer peaks com- 
pletely disappear!) 

Thus, we resort to a general alternative method to esti- 
mate i?coox(r), which requires to scan x{H), as defined in 
Eq. (|34|). as a function of the field H . This method would 
have been very inconvenient at low temperatures - the 
increments 5H of the field H necessary to perform such 
a scan would need to satisfy the condition 5H <C AH 
and since AH is so small {Eq. (p6| } and HcoexiT) is 
not known beforehand, an enormous (and not reason- 
able) effort of computer resources would be implied (and 
furthermore the transition would easily be missed due to 
hysteresis). However, for T > Tf){L,D) hysteresis is no 
longer a severe problem, and the rounding of the tran- 
sition is much smaller, since now (Nc is the number of 
spins in a cross sectional plane, N = NcL) 



Xmax(r) « mD{T)N, 



AT) 



(37) 



and thus Xmax{T) is much smaller than in the region 
where Eq. ([Ml) holds. So AH is no longer so small; fur- 
thermore, one can get a first estimate for iJcoex(T') by an 
extrapolation from the region where Pl^d{M) has only 
two peaks, and the analysis as described above works. 

Thus, we have scanned the region of interest choos- 
ing steps AH = 0.0002, carrying out runs with 2 million 
MC steps per spin at each state point {H, T). Histogram 




FIG. 14. Plot of x' vs. H for the case ffi = 0.75, i? = 4, L = 
300 and increasing temperatures from left to right obtained by 
explicit simulations at each value of H . b) Plot of x' vs. H for 
the case Hi = 0.75, R — 5, L = 200 and several temperatures 
as indicated using extrapolation to different values of H. 



reweighting methods were used to improve the accuracy 
of the results, as is standard^. Fig. [T4k presents typical 
data for the case R = A, L = 300. Since we have found 
that also x' ^ts defined in Eq. pip has a sharp peak at 
H « Hcoox we used the location Xmax ^ well. Note that 
in the region where P]^ jj(M) has a single peak, the width 
of this peak is rather narrow if H differs appreciably from 
HcoexiT), since there then the state of the pore is uni- 
form, no nucleation of two-phase fluctuations takes place. 
Then (M^) sa (l-M^I)^ irrespective of the value of the peak, 
and x' is of order unity. Only for H near -ffcoox will the 
distribution Pl,d{M) show some anomalous broadening, 
resulting from the fluctuations associated with the co- 
existence of multiple domains. Therefore, recording the 
maxima of x'j which at temperatures near Tq(L,D) are 
much easier to sample (Fig. [Hb). is a useful method to 
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FIG. 15. Configuration snapshots of a system with Hi = 
0.75, 7? = 4, L = 200 at Hcoc^T) for T = 2.5 (a) and T = 
3.65 (b). Each snapshot shows for each peak of Pl,d{M) an 
outside projection (upper part) and a cross section containing 
the cylinder axis (lower part). Negative spins are shown in 
black, positive spins in yellow. Note that for T = 2.5 one has 
two peaks for Pl,d{M,T), while for T = 3.65 one has three 
peaks. The snapshots for the middle peak in part (b) show 
the multiple domain structure clearly. 



estimate i?coex(T) for T > To{L,D), see Fig. [Hj. 

It is interesting to note that the coexisting phases of 
the cylinder (in the region T < Tq(L,D)) are inhonio- 
geneous. This is evidenced both by snapshot pictures 
(Fig. [T5]) and radial order parameter profiles (Fig. [TB)) 
taken for our systems. One can see that on the outside 
surface of the cylinder (seen in the projection snapshots 
of Fig. lisp there is always more disorder. The reduction 
of the local magnetization at the surface, when the bulk 
of the cylinder has positive magnetization, can be inter- 
preted as a precursor of wetting phenomena. Of course, 
true wetting layers cannot form in nanopore cylinders, 
and hence we also do not find a transition as proposed by 
Liu et ali^. Crossing the wetting transition temperature 
T^, (_ffi it was predicted that a transition from 

"plugs" to "capsules" should occur^*", and an attempt 
was madeii^ to locate this transition by Monte Carlo 
simulations in L± x L± x L systems with L± = 14, 20 
and 28, varying L from L = 40 to L = 320, and various 
choices of Hi. However, for such rather wide and short 
pores the problem of multiple domains did not yet arise, 
and the issues about intrinsic rounding of transitions in 
the quasi-one-dimensional pore geometry were not stud- 
ied in these investigation s^^' 

As a final point of this section, we present in Fig.lTTlthe 
"phase diagram" of our model, both in the T — M plane 
and the H — T plane. Note that this "phase diagram" is 
only meant to describe the phase coexistence that per- 
sists if £,d{T) D on a local scale in long cylinders (or 
in short cylinders, if £,d{T) still exceeds L). On the scale 
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FIG. 16. a) Magnetization profiles m(r) as a function of the 
distance r from the cylinder axis, for Hi = 0.75, R = 4, 
L = 200, and T = 1.5, 2.5 and 3.5, as indicated. The upper 
set of curves belongs to states with positive magnetization 
of the pore, the lower set for states with negative magneti- 
zation, for H = -H'coGx(r). In all cases Pl,d{M) still has 
only two peaks, b) Magnetization profiles m(r) as a function 
of the distance r for T — 3.65, showing both profiles from 
the peaks representing the "pure phases" and from the right 
and left parts of the middle peak. These latter profiles were 
extracted from local slices through the cylinder (the contri- 
bution of slices containing the interfaces turn out to be still 
negligible at this temperature.) 



of the axes chosen in Fig. [T71 the "phase boundaries" 
still look sharp, although the transition line in Fig. [TTb 
is intrinsically rounded over some width A.H, but for the 
temperature region shown, the rounding is still small. 
However, this phase diagram cannot uniquely be contin- 
ued up to a "capillary critical point" : when (T) has 
decreased to a value comparable to D, the rounding gets 
very strong, and different criteria to locate a "transition" 
will no longer coincide (e.g., the position of a maximum 
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FIG. 17. "Phase diagram" of the cyhndrical Ising pore for the 
case R — 4, Hi — 0.75, and three choices of L plotted in the 
(T,M) plane (a) and in the (H,T) plane (b). In part (a), data 
for Mcoex(T) for the spontaneous magnetization of a three 
dimensional bulk Ising model are included. The data for the 
two coexisting phases in (short) cylinders are extracted from 
the positions of the two peaks at H — //cocx(r), respectively. 



for x{H) and x'{H) will no longer agree, etc.) Thus, 
our "phase diagram" ends in an "open way" at T = 4.0: 
for pores as narrow as i? = 4, the difference between 
the order parameter of the two coexisting phases then 
has already decreased significantly, and at slightly higher 
temperatures it is no longer possible to distinguish the 
regions of the "pure coexisting phases" inside the pore 
from the interfaces separating them. However, it is al- 
ways of interest to study in very long pores the transi- 
tion at Cd(T) ~ L/3 from the multiple domain phase 
coexistence at Hcocx{T) to the "pure" coexisting phases 
inside the pore. Fig. [TSk compares this transition for 
three choices of Hi: we see that increasing the surface 
field shifts the transition to lower temperatures, but the 



FIG. 18. a) Weight W of the central peak oi Pl,d{M) plotted 
vs. temperature for Ising pores of radius _R = 4 and three 
choices of the surface field, Hi = 1.5, Hi = 0.75 and Hi = 
(from left to the right in different colors). In each case the 
pore lengths L = 600, 400, and 200 are included (from the 
left to the right with different line styles), b) Magnetization 
of Ising cylinders of radius R — 4, length L — 200, surface 
magnetic field Hi — 0.75, plotted vs. field H — Hcoex{T) at 5 
temperatures, as indicated. Bulk field H was varied in steps 
of 0.001. 



qualitative characteristics stay the same. 

As in the case of the Ising strip, we can verify that 
in the same region of temperatures where the change of 
Pl.d{M) from the double-peak distribution to the triple- 
peak distribution occurs (cf. Fig. [T3k ) the hysteresis in 
the magnetization process vanishes (Fig. 118b ). namely 
near T sa 3.6. The figure shows that at this temperature 
one still can identify clearly the difference in order param- 
eter of the vapor-like and liquid-like branch of the lattice 
gas model. The data in Fig. [T5b ) are for a rather short 
pore, namely L = 200: It is clear (cf. also Fig. [T5b ) that 
for a longer L this change of Pl.d{M) occurs for lower 
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temperatures, and also the onset of hysteresis occurs at 
the lower temperature the larger the pore length L is 
considered. 



III. COLLOID-POLYMER MIXTURES CONFINED IN 
CYLINDERS: A MONTE CARLO STUDY OF THE 
ASAKURA-OOSAWA MODEL 

Colloidal dispersions have become model systems for 
the study of phase behavior of condensed matter, since 
the large size of the colloidal particles allows the use 
of experimental observation techniques that cannot be 
used for small molecular systems. Furthermore, inter- 
actions among colloidal particles are tunable to a large 
extentii^^-i^. Colloid-polymer mixturesi22r-i21 have been 
particularly suitable to study liquid- vapor-like phase sep- 
aration into colloid-rich and colloid-poor phases, includ- 
ing their interfacial behavior. There also exists a very 
simple theoretical model, due to Asakura and Oosawai^S 
and Vriji^ (henceforth referred to as "AO model" ) , well 
suited for Monte Carlo simulation studie o^'^'^^^" -^. In 
this model, colloids are simply described as hard spheres 
of radius Rc while polymers are soft spheres of radius Rp. 
While overlap among colloids and between polymers and 
colloids is strictly forbidden, i.e. the potential energy is 
given by 



Ucc{r < 2R,) = C50, Ucc{r > 2R, = 0) 



(38) 
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Upc{r <Rp+Rc)^^, Upc{r>Rp + Rc} = , (39) 

two polymer coils can interpenetrate and hence overlap 
with no energy cost, Upp{r) — irrespective of distance 
r. Vink et alji^"-i22 have already performed an exten- 
sive study of capillary condensation for colloid-polymer 
mixtures confined between two parallel hard walls a dis- 
tance D apart, and have shown that for distance D of 
the order of a few colloid diameters a crossover from 
three-dimensional to two-dimensional Ising critical be- 
havior occurs. Due to the large sizes of colloid particles, 
it should be experimentally feasible to also study capil- 
laries which are only a few colloids' diameters wide, and 
since for particles in the size range of a /^m the atomistic 
corrugation of real walls clearly is negligible, fairly ideal 
conditions should be realizable. 

In the present work, we have extended this worki22r-i^ 
to confinement in cylindrical pores of diameters D — 
12i?c and lengths L up to L = 540i?c, for Rp/Rc = 0.8. 
In the following, Rc — 1 shall be used as unit of length in 
this section. The Monte Carlo simulations were carried 
out in the grand-canonical ensemble, choosing the chem- 
ical potential of the colloids and the polymer reservoir 
packing fraction rf^ 

r,; = {4TTRl/3)eMPp/kBT) , (40) 



FIG. 19. Radial density of colloids {pc{r)) and polymers 
{pp{r)) plotted vs. distance from the cylinder axis for r^p = 
1.30, L — 60, D = 12. Case a) shows these profiles in the 
vapor-like phase, case (b) for the liquid-like phase. 
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FIG. 20. Parallelization scheme of the "Wang-Landau" algo- 
rithm. The first row of numbers is the CPU index. Mi is the 
number of Monte Carlo steps, which is performed. Wi is the 
weight function of GPU i. The brackets denote an average 
weighted by the number of MC steps. The average replaces 
the weight function Wi of every process i. 
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where /Xp is the chemical potential of the polymers, 
as independent control variables. Observables of inter- 
est then are both global average densities pp — Np/V, 
Pc = Nc/V of colloids and polymers {Np, Nc are the 
total number of polymers and colloids in the volume 
V = [ttD'^ /A)L as well as the corresponding radial den- 
sity profiles Pp{r), pdr), see Fig. [121 One can see that in 
the vapor-like phase the polymer density is reduced near 
the pore wall, while colloids are attracted to the pore 
wall both in the vapor-like and liquid-like phase. Note 
that phase coexistence in the pore was located as for the 
thin film case by scanning the chemical potential pc un- 
til one finds a double-peak distribution, where then the 
equal weight rul o^^^i^^'^ is applied to estimate the value 
of the chemical potential at coexistence /ic, coox- As for 
the case of the AO model in the bulki2iii2^ and in thin 
film geometryi^""-i^, cluster moves^^ and successive um- 
brella sampling methods-^- are applied throughout. For 
large systems a parallel version of the "Wang-Landau" 
algorithmic was implemented. The idea, schematically 
shown in Fig. [501 is to correlate a priori independent 
simulations by taking averages over the weight functions 
iteratively generated by the "Wang-Landau" algorithm 
on each CPU. The average between the weight functions 
of the single simulations is weighted with the MC steps 
done so far and is used as a new weight function for all 
CPUs. This procedure uses only a very small amount of 
communication and scales therefore almost linearly up to 
1440 CPUs. In comparison to the same number of non- 
communicating independent simulations, a reduction of 
the systematic error of the biasing algorithm was ob- 
served, which leads to a faster convergence of the weight 
function to the true free energy landscape. 

While for the case shown in Fig. [THl the state of the 
cylinder at phase coexistence is axially homogeneous, 
and this fact also shows up in the probability distribu- 
tion Pl,d{i1c), Vc — {4:TrRl/3)pc being the colloid packing 
fraction. Fig. [5T1 since Pl,d{i1c) just has two peaks and 
is flat in between, at lower values of ijp one again finds a 
distribution with three peaks. As in the Ising case, the 
interpretation of the distribution exhibiting a "central" 
peak is the formation of multiple domain walls across the 
pore (Fig. [221). 

Fig. [55b shows the average colloid density {pc) as a 
function of /Xc- The maximum value of the fluctuation 
{{pc— {Pc)Y) studied as a function of rfp for several choices 
of L is shown in Fig. [23b . Also the corresponding phase 
diagram is shown (Fig. [55b). While in the bulk well- 
defined vapor-liquid like phase coexistence occurs, ending 
in a critical point at rfp ^r — 0.765, ?7c,cr = 0.13, the 
phase coexistence in the cylindrical pore exists over a 
finite correlation length S^o only. The value of {rjc) in the 
coexisting vapor-like and liquid-like phases depend on L 
only very weakly: however, the larger L the larger rfp has 
to be chosen to ensure that one still has phase coexistence 
between single-domain states in the pore, rather than a 
multiple domain structure. 

Again it is of central importance to verify the connec- 
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FIG. 21. Probability distribution Pl,d(t]c) of the colloid pack- 
ing fraction in a cylinder of diameter D = 12 and length 
L = 180 (a) and for _D = 6, L = 100 (b). 



FIG. 22. Snapshots for the system with D = 12 and L = 540 
at the polymer reservoir packing fraction 77^ = 1.15. Colloids 
are shown in yellow. The lower picture shows the cylindrical 
simulation box sliced at the plane a; = 0, while the upper pic- 
ture visualizes the projection of the particles at the confining 
border. 



tion between the change of the distribution P{rjc) with 
decreasing rf^ from the double peak behavior at large 
rfp to the three-peak behavior at somewhat smaller rf^ 
(cf. Fig. [21]) and the disappearance of hysteresis at a value 
of rfp which is still distinctly larger than the pore criti- 
cal temperature (where in Fig. 123b the vapor-like and 
liquid-like branches of the coexistence curve of the fluid 
confined in the pore have merged). This connection is 
verified by Fig. [Ml which shows that for D = 12 and 
L = 180 hysteresis indeed disappears in between = 1.2 
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FIG. 23. a) Average density (pc) of the colloids in a pore of 
linear dimensions L — 60, D — 12 plotted vs. fic at various 
values of r/p as indicated, b) Plot of the maximum value of 
the density fluctuation for D — 12 and various L as indicated, 
plotted vs. rjp. c) Phase diagram of the AO model in the plane 
of variables {rip, rjc) shown for cylinders of diameter D = 
12 and three choices of L. Full curve shows the coexistence 
curve for the corresponding bulk AO model. The symbols at 
Tjc ~ 0.16 to 0.17 show the transition from the single-domain 
to the multiple domain state in the pore, d) Barrier AF/T 
against nucleation of interfaces in the AO model confined to 
a cylindrical pore of diameter _D = 12 plotted versus inverse 
polymer reservoir packing fraction. 
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FIG. 24. Two hysteresis plots for the AO model. The chemi- 
cal potential was varied in steps of O.OOlfcsT. Several simula- 
tion runs (up to 38) were averaged. For high polymer reservoir 
packing fractions large sample to sample fluctuations occur. 
The open symbols show data for which the chemical poten- 
tial was increased step-wise while the full symbols show data 
for which the value of the chemical potential was decreased 
step by step, (a) shows the disappearance of the hysteresis 
for a system with D — 12 and L = 180. (b) shows the dis- 
appearance of the hysteresis for a system with D = 6 and 
L = 100. 



and rjp = 1.1, while the eoexistence curve branches exist 
up to about rjp = 1.0 (Fig. [^5b). Again we predict that 
this difference between the pore critical temperature and 
the hysteresis critical temperature should increase with 
L. As an analogue to Fig. [8] the free energy barriers are 
shown in Fig. [23J1 for various choices of L. 

An alternative way to explore this transition from ax- 
ially symmetric phase coexistence in the cylindrical pore 
to a multiple domain structure uses a very long pore 
{L = 1800) which is cut into a one-dimensional array of 
Ns subsystems, and recording the distribution Pd.n^ i^c) 
of the number of colloids in the subsystems (Fig.BSj). One 
can nicely see that for short enough subsystems (i.e., for 
Ng > 60) the subsystem is typically homogeneous, since 
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FIG. 25. Three-dimensional plot (left part) and contour plot 
(right part) of Pd,Ns{Nc) vs. number of subsystems Ns, for 
rjp = 1.1 (first two graphs) and rip = 1.5 (last two graphs). 
Note that the number of colloids Nc is normalized by the 
length L/Ns. 



there occur just two peaks with a minimum in between. 
However, for very large L/Ng one still finds the mid- 
dle peak, as a signature of the multiple domain struc- 
ture, and the transition between both types of behaviors 
(as a function of L/Ng or 77^, respectively) is completely 
gradual. Thus, the fact that the coexisting phases in 
the phase diagram of Fig. [23b show practically no L- 
dependence, and the fact that the equilibrium isotherms 
(Fig. [23k ) at large 77^ have an almost perpendicular part 
should not be taken as evidence that in the cylindrical 
pore a sharp, well-defined phase transition exists: as in 
the Ising model, the transition is rounded, but for large 
the extent of rounding is very small. 



IV. CONCLUSIONS 

In this paper, the characteristic features of phase tran- 
sitions of Ising-like systems in a quasi-one dimensional 
geometry have been explored by Monte Carlo simula- 
tions for four generic models: (i) Ising L x D strips with 
D and periodic boundary conditions throughout (ii) 
Ising "cylinders" of length L and cross section containing 
Ncr sites enclosed by a circle of radius i?, with a "missing 
neighbor" boundary condition that does not destroy the 
symmetry between the coexisting phases in the ground 
state; (iii) the same model as in (ii), but with a surface 
field Hi acting on the spins which have "missing neigh- 
bors" , so that the spin reversal symmetry is broken, and 
the model (interpreted as a lattice gas) exhibits capil- 
lary condensation; and (iv) the AO model confined to 
cylindrical pores of diameter D and length L, confined 
by hard repulsive walls, as an off-lattice model that lacks 
any particular symmetries already in the bulk. 

We have shown that all models exhibit qualitatively 
similar behavior, namely two strongly rounded tran- 
sitions occur when at phase coexistence conditions 
the temperature (or temperature-like variable, such 
as {np)~^ in the case of the colloid-polymer mixture, 
respectively) is lowered: at a temperature rather close to 
the critical temperature of the bulk, a rounded transition 
occurs from the disordered phase (which is axially sym- 
metric but may have nontrivial order parameter profiles 
in the phase perpendicular to the cylinder axis, induced 
by the boundaries, if there is no complete symmetry 
between the coexisting phases) to a locally ordered 
phase, where the size of the domains ^.oiT) in axial 
direction exceeds distinctly the pore diameter, so that a 
long cylinder (L ^ (,d{T) ^ D) is characterized by a 
sequence of interfaces across the cylinder axis. The order 
parameter distribution at coexistence is then a very 
broad Gaussian characterized by a very large response 
function (if the transition is studied as a function of the 
field conjugate to the order parameter, the rounding 
is exponentially small in the cross-sectional areas of 
the cylinder). At i/3 w ^d{T), i.e. at T = To{L,D), 
a second, again rounded, transition occurs, where the 
state of the system is again axially uniform and either 
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the vapor-like or liquid-like phase dominates. When one 
studies the kinetics of the transition between vapor-like 
and fluid-like phases, varying the field conjugate to 
the order parameter, one finds pronounced hysteresis 
in this low temperature region, T < To{L,D) where 
^£)(T) > L, but these hysteresis loops get narrow when 
T « To{L,D) and vanish completely for T > To{L,D). 
Thus, we suggest that the "hysteresis critical point" 
can be associated with the lower temperature Tq{L,D) 
rather than the upper pseudo-critical temperature of the 
capillary (where ^d{T) ~ D and the difference in order 
parameter between the coexisting phases disappears). 
A prediction that could be tested experimentally is 
our result that the difference between this "hysteresis 
critical temperature" and the "pore critical tempera- 
ture" should increase logarithmically with the length 
L of the cylindrical pore. We hope that our study 
stimulates additional experimental work using pores of 
both well-controlled diameter and length to check our 
predictions and thus confirm that a long-standing puzzle 
about the phase behavior of fluids adsorbed in pores is 
now better understood. 
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